



//The first thing I do is calculate the total MW reductions from the policy when extended

//Step 1 - figure out customer counts from demographic data

use "stata data/TOU_demographic 2", clear
drop if stopdate < mdy(6,1,2014)
drop if sastart >= mdy(6,1,2014)



gen class = ""
replace class = "A1" if inlist(rs,"A1","A1F","A1X","HA1","HA1X")


foreach type in E19P E19S E19V E20P E20S E20T A6 A10{
gen `type' = regexm(rsched,"`type'")
replace class = "`type'" if `type' == 1
}

replace rs = "" if class !=""

drop if rs != ""

duplicates drop spid, force

rename climate climate_str
encode climate_str, gen(climate)

gen coastal = 0 
replace coastal = 1 if inlist(climate,5,6,2)

preserve
gen counter = 1
collapse (sum) counter, by(class coastal)
rename class rs
save "stata data/Count of customer class in PGE terms by coastal", replace
restore


gen counter = 1
collapse (sum) counter, by(class)
rename class rs
save "stata data/Count of customer class in PGE terms", replace




//Step 2 - figure out how many inland customers have between 4000 and 50000 consumpitoon
/// Percent of customers in ranges
use "stata data/Sa rs for 32 to 45 april 11 2016", clear

merge m:1 sa using "stata data\TOU_demographic 2 climate"
keep if _merge == 3
drop _merge

gen coastal = 0
replace coastal = 1 if inlist(climate,5,6,2)

keep if coastal == 0

keep if inrange(bin,-27,27)
merge 1:1 sa using "stata data/larger total_kwh file"
keep if _merge == 3

count if inrange(total,4000,50000)
local in_range = `r(N)'
count
local total = `r(N)'

gen percent_in_range = `in_range' / `total'
 

keep percent_in_range
duplicates drop
gen rs = "A1"
save "stata data/percent of customers in 4000_50000 range", replace


//step 3 - calculate opt out rate for people placed in for first wave




use "stata data/extract file september 2016", clear

gen date_str = substr(pdp_start_date,1,10)
gen pdp_start = date(date_str,"YMD")

gen end_date_str = substr(pdp_end_date,1,10)
gen pdp_end = date(end_date_str,"YMD")

drop end_date_str date_str

format pdp_start pdp_end %td

gen pdp_start_year = year(pdp_start)

gen optionally_enrolled = 0
replace optionally_enrolled = 1 if enroll_optional !=""
drop enroll_optional 

rename sa_id sa

merge 1:1 sa using "stata data/TOU_demographic 2"
keep if _merge == 3

keep sa  status_name action_type action_status program_status optionally_enrolled pdp_start* pdp_end rs stopdate

gen pdp = 0
replace pdp = 1 if program_st == "Enrolled"

//I am defining defaulted as placed on to opt out rate. 
//defaulted is anyone that has an action type. People who opt out are inactive status 

gen opt_out = 0
replace opt_out = 1 if action_type == "Optout"
replace opt_out = 1 if action_type == "Unenroll"

keep if rs == "HA1X"
keep if pdp_start_year == 2014
drop if optionally_enrolled == 1

save "temp/opt out new extract updated", replace


use "temp/opt out new extract updated", clear

drop if action_status == "Closed SA" //do this so I"m not counting the people that just seem to drop out


gen dropout_date_ym = mofd(pdp_end)


gen dropout_before_first_year = 0
replace dropout_before_first_year = 1 if dropout_date_ym < 665 & dropout_date_ym !=.


gen dropout_during_first_summer = 0
replace dropout_during_first_summer  = 1 if ///
	inrange(dropout_date_ym,665,669) & dropout_date_ym !=.

gen dropout_after_first_summer = 0
replace dropout_after_first_summer  = 1 if ///
	dropout_date_ym>=670 & dropout_date_ym !=.

gen dropout_first_year = dropout_before_first_year
replace dropout_first_year  = 1 if dropout_during_first_summer == 1
	
collapse (mean) dropout_before dropout_during dropout_after opt_out dropout_first_year



rename opt_out opt_out_ratio	
gen rs = "A1"

save "stata data/opt out ratio", replace



//step 4 - calculate total MW savings
//customer count for inland firms
use "stata data/Count of customer class in PGE terms by coastal", clear
keep if rs == "A1" & coastal == 0
save "stata data/inland count", replace


//pull in temp bin regression results
use "stata data/regsave/FE_new_all_tempbins_inland_event", clear
drop if inlist(var,"temp_f","temp_f_2")
drop if pval >.05
keep var coef
split var, p(_bin_)
destring var2, gen(bin)
keep bin coef
gen id = 1
reshape wide coef, i(id) j(bin)
rename coef3 coef_100_101
rename coef4 coef_102_104
drop id
gen rs = "A1" //what I merge with
save "stata data/coef FE with temp bins", replace


//pull in non-event days
use "stata data/regsave/FE_new_all_inland_non_event", clear
keep if var == "pdp"
gen rs = "A1"
keep coef rs
rename coef coef_non_event
save "stata data/coef_non_event", replace



use "stata data/regsave/FE_new_all_inland_event", clear
gen rs = "A1"

keep if var == "pdp"



merge 1:1 rs using "stata data/coef FE with temp bins"
drop _merge

merge 1:1 rs using "stata data/coef_non_event"


gen percent_reduction_iv =  -( 1-exp(coef))  
gen percent_reduction_non_event =  -(1-exp(coef_non_event))

gen firm_level_reductions = avg_kwh * percent_reduction_iv
gen firm_level_reductions_non_event = avg_kwh * percent_reduction_non_event


// make reductions for the different temp bins. Anything below with 100 or 102 on it is the same new addition. 
gen percent_reduction_iv_100_101 =  -( 1-exp(coef_100_101))
gen percent_reduction_iv_102_104 =  -( 1-exp(coef_102_104))
gen firm_level_reductions_100_101 = avg_kwh * percent_reduction_iv_100_101 
gen firm_level_reductions_102_104 = avg_kwh * percent_reduction_iv_102_104


keep coef avg_kwh firm_level_reductions* rs percent_reduction*


merge 1:1 rs using "stata data/inland count"
drop coastal _merge


merge 1:1 rs using  "stata data/percent of customers in 4000_50000 range"
drop _merge

merge 1:1 rs using "stata data/opt out ratio"

gen customers = counter * percent_in_range

gen total_reductions = firm_level_reductions * customers *(1-opt_out_ratio)
gen total_reductions_non_event = firm_level_reductions_non_event * customers*(1-opt_out_ratio)


//do temperature specific savings 
gen total_reductions_100_101 = firm_level_reductions_100_101 * customers *(1-opt_out_ratio)
gen total_reductions_102_104 = firm_level_reductions_102_104 * customers *(1-opt_out_ratio)



keep total_reductions*	 avg_kwh percent_reduction* firm_level_reductions* customers opt_out_ratio

gen year = 2015
gen mw_reductions_fe = total_reductions / 1000
gen mw_reductions_non_event = total_reductions_non_event / 1000
gen mw_reductions_fe_100_101 = total_reductions_100_101 / 1000
gen mw_reductions_fe_102_104 = total_reductions_102_104 / 1000

keep year mw* avg_kwh percent_reduction* firm_level_reductions* customers opt_out_ratio

save "stata data/total reductions", replace


///this next step is a bit odd. I previously was using this data as part of my calculations, but now I'm just using the date structure 

foreach month in june july august september october {
insheet using "public raw data\\`month' 2015.csv", comma clear
save "stata data//`month' 2015 demand raw", replace
}

use "stata data//june 2015 demand raw", clear
foreach month in july august september october {
append using "stata data//`month' 2015 demand raw"
}
 
 
keep if market == "RTM"
gen hour = opr_hr - 1
keep if label =="RTM 5Min Load Forecast"
keep if  tac_area_name == "PGE-TAC" //now cut it down to just PGE


gen minute = (opr_interval - 1)*5
order opr_interval minute



gen date = date(opr_dt,"YMD")
keep hour mw date opr_interval  minute tac_area
format date %td

gen month = month(date)
gen year = year(date)
gen day = day(date)

gen double date_clock = mdyhms(month,day,year,hour, minute,0)
format date_clock %tc


//this is a strange thing. Apparently there are just some missing observations in the data for a few 5 minute increments
//so I just used the thing before since it had to be close
//below is the proess I do that by. I just expand out and add in 300000 milliseconds, or 5 minutes. Seems to work well
sort date_clock //this is important to have in there!
gen time_difference = date_clock - date_clock[_n+1]
gen expander = -time_difference/300000
expand expander

sort date_clock

bysort date_clock: gen within_clock_counter = _n

gen double date_clock_adjusted = date_clock + (within_clock_counter -1)*300000
format date_clock_adjusted %tc
order date_clock*

drop date_clock
rename date_clock_adjusted date_clock

//this bit is b/c when I do the adjustment, sometimes if you have a 13 hour spread into a 14 hour, the 
//hour won't update. So I just repull it from the date clock
drop hour
gen hour = hh(date_clock)

save "stata data/5 minute demands by region before cut", replace

use "stata data/5 minute demands by region before cut", clear




merge m:1 date using "stata data/2015 pdp event dates"
drop if _merge == 2
gen event_day = 0
replace event_day = 1 if _merge == 3
drop _merge

gen dow = dow(date)

gsort -mw



///grab the MW reductions from up above
//year merge is just to get it onto every observation
merge m:1 year using "stata data/total reductions"
drop _merge


save "stata data/side for super-peak figure", replace





foreach p_val in  1.1 1.5 .6 {
foreach mw_scenario in  "mw_reductions_fe_100_101" "mw_reductions_fe_102_104" "mw_reductions_fe" {
*** BELOW A LOOP OVER DIFFERNT SCENARIOS 
use "stata data/side for super-peak figure", clear


global new_price_change_varlabel = `p_val'*10
global new_price_change `p_val'


global mw_reductions_value `mw_scenario' //I set it up as a global so I can name the file...



gen Q_1 = avg_kwh + firm_level_reductions

*replace Q_1 = avg_kwh + firm_level_reductions_102_104  //< MANUAL CHANGE -- keep this off unless doing a 102_104 robustness 

gen epsilon_ces = ln(avg_kwh/Q_1)/ln(.24756/.84756)

*now derive A at the anchor point of the initial starting position
*A = Q_0/p_0^epsilon

gen A = avg_kwh/(.24756^epsilon_ces)


global sensitivity_scalar 1 // This is no longer used in main paper but I'm leaving in  case sensitivity is needed

*quick CES break
gen firm_reductions_ces = avg_kwh - A*(.24756+$new_price_change)^epsilon_ces

gen mw_reductions_ces =  (firm_reductions_ces  * customers*(1-opt_out_ratio))*$sensitivity_scalar/1000


//this is new code added 6/3/21 - it allows me to swap in mw_reductions_101_102 in for mw_reductions, which is the default
//when I sawp in the new value, that's how I can see effects with higher reductions
//can swap in mw_reductions_100_101 or mw_reductions_102_104s

gen mw_reductions = -1* $mw_reductions_value *$sensitivity_scalar //sensitivity scalar is usually 1, can make 2 for robustness check


//here I can put in alternate savings. 
//I calculate the SLOPE as -.6/156 (price change over mwh change) is around .00384
//so I can solve for a new quantity as the new price change divided by that slope


gen new_price_change = $new_price_change //this is where I can choose a new price change. Set to $1 for now

gen demand_slope = .6/mw_reductions
gen demand_slope_sim = .6/mw_reductions



gen mw_reductions_sim = new_price_change/demand_slope_sim
gen mw_reductions_sim_ces = (avg_kwh - A*(.24756+$new_price_change)^epsilon_ces) * customers*(1-opt_out_ratio)*$sensitivity_scalar/1000




//here I want to derive how much better off people are from the 1 cent reduction in all other hours
gen mw_increase_1cent = .00977/demand_slope //<- increase in MW on grid from decrese in price in all other hours
gen mw_increase_1cent_sim = .00977/demand_slope_sim //<- increase in MW on grid from decrese in price in all other hours


save "stata data/in progress 9485", replace


use "stata data/in progress 9485", clear


gen welfare_losses  = mw_reductions  * $new_price_change  *.5 *1000 / 12  //NOTE: the *1000 is to go from MW to KW, and the 12 is to go from hours to 5 minute chunks
replace welfare_losses  = 0 if hour < 14
replace welfare_losses  = 0 if hour >17


gen welfare_losses_sim  = mw_reductions_sim * $new_price_change *.5 *1000/12
replace welfare_losses_sim  = 0 if hour < 14
replace welfare_losses_sim  = 0 if hour >17

gen Q_sim = A*($new_price_change + .24756)^epsilon

gen welfare_losses_R1_tall_triangle = mw_reductions  * 1.1  *.5 *1000 / 12  //on 6/7/21  R1 asked to see the MW reductions from .85 but the welfare losses at the 135 price. it requires recalculating this new tall triangle with the 
///1.1 addeder (i.e. 1.35) price but same 118 mW reduction. This will only work for default scenaro
replace welfare_losses_R1_tall_triangle  = 0 if hour < 14
replace welfare_losses_R1_tall_triangle  = 0 if hour >17


///REMEMBER this functional form is flipped on the axes. So you subtract the area to the left of the curve. Which is why it's .6*Q1 you are subtracting
gen welfare_tri_per_person_ces = (A/(epsilon_ces +1))*((.24756+.6)^(epsilon_ces+1)-(.24756)^(epsilon_ces+1))-(.6)*(Q_1)
gen welfare_tri_per_person_sim_ces = (A/(epsilon_ces +1))*((.24756+$new_price_change)^(epsilon_ces+1)-(.24756)^(epsilon_ces+1))-($new_price_change)*(Q_sim)


//so the customers*(1-opt_out_ratio)/1000 part converts it up to MW and the *1000 converts it back to kwh. I do this (I know it's confusing)
*to better match what's going on here. 
gen welfare_loss_triangle_ces = welfare_tri_per_person_ces * customers*(1-opt_out_ratio)/1000 *1000 / 12
gen welfare_loss_triangle_sim_ces = welfare_tri_per_person_sim_ces * customers*(1-opt_out_ratio)/1000 *1000 / 12

//these are the same, just used different notation to make it clearer above
gen welfare_losses_ces = welfare_loss_triangle_ces 
replace welfare_losses_ces  = 0 if hour < 14
replace welfare_losses_ces  = 0 if hour >17


gen welfare_losses_sim_ces = welfare_loss_triangle_sim_ces
replace welfare_losses_sim_ces  = 0 if hour < 14
replace welfare_losses_sim_ces  = 0 if hour >17





global MC_peaker .10266 //<-----this is where I enter the marginal cost of energy from peakers

//quickly add in half hour stuff
gen half_hour_adder = 0
replace half_hour_adder = 1 if inrange(minute,0,25)
gen half_hour = hour * 2
replace half_hour = half_hour + half_hour_adder

gen period = "off peak"
replace period = "part peak" if inrange(half_hour,17,23)
replace period = "part peak" if inrange(half_hour,36,42)
replace period = "peak" if inrange(half_hour,24,35)

gen tou_price = .
replace tou_price = .24756 if period == "peak"
replace tou_price = .23875 if period == "part peak"
replace tou_price = .21195 if period == "off peak"


save "stata data/progress 999888", replace


use "stata data/progress 999888", clear


///
///Now do the 1 cent decrease. Rememeber there is no sim value here b/c I assume it's always a 1 cent refund. 
///

//linear
gen welfare_triangle_gains_1cent = mw_increase_1cent * .00977 * .5 *1000 / 12  //NOTE: the *1000 is to go from MW to KW, and the 12 is to go from hours to 5 minute chunks
gen welfare_rect_gains_1cent = (tou_price-.00977 - $MC_peaker) * mw_increase_1cent * 1000/12

gen welfare_triangle_gains_1cent_sim = mw_increase_1cent_sim * .00977 * .5 *1000 / 12  //NOTE: the *1000 is to go from MW to KW, and the 12 is to go from hours to 5 minute chunks
gen welfare_rect_gains_1cent_sim = (tou_price-.00977 - $MC_peaker) * mw_increase_1cent_sim * 1000/12



//CES
gen mw_increase_1cent_ces = -1*(avg_kwh - A*(.247560-.00977)^epsilon_ces)* customers*(1-opt_out_ratio)/1000
label variable mw_increase_1cent_ces  "Increase in consumption in off hours due to 1 cent decrease - CES"

gen welfare_1cent_per_person_ces = ((A/(epsilon_ces +1))*((.24756)^(epsilon_ces+1)-(.24756-.00977)^(epsilon_ces+1))-(.00977)*(avg_kwh))

gen welfare_triangle_gains_1cent_ces  = welfare_1cent_per_person_ces *customers*(1-opt_out_ratio)/1000 *1000 / 12
gen welfare_rect_gains_1cent_ces = (tou_price-.00977 - $MC_peaker) * mw_increase_1cent_ces * 1000/12



//add up welfare gains from 1 cent increase
gen welfare_gains_1cent = welfare_triangle_gains_1cent + welfare_rect_gains_1cent 
replace welfare_gains_1cent  = 0 if inrange(hour,14,17) & event_day == 1

gen welfare_gains_1cent_sim = welfare_triangle_gains_1cent_sim + welfare_rect_gains_1cent_sim
replace welfare_gains_1cent_sim  = 0 if inrange(hour,14,17) & event_day == 1



gen welfare_gains_1cent_ces = welfare_triangle_gains_1cent_ces + welfare_rect_gains_1cent_ces
replace welfare_gains_1cent_ces  = 0 if inrange(hour,14,17) & event_day == 1

save "stata data/move on to rectangle below triangle process", replace


use "stata data/move on to rectangle below triangle process", clear

///
///Now move back to welfare losses of the rectangle below the main triangle
///


//these take into account the difference between the MC of generation and the retail price charged
gen welfare_loss_rectangle = (.24756 - $MC_peaker) * mw_reductions * 1000/12
gen welfare_loss_rectangle_sim = (.24756 - $MC_peaker) * mw_reductions_sim * 1000/12

gen welfare_loss_rectangle_non_event = (.24756 - $MC_peaker -.01) * -mw_reductions_non_event * 1000/12 //I added the -.01 b/c there was a discount already...
replace welfare_loss_rectangle_non_event  = 0 if event_day == 1

gen welfare_loss_rectangle_ces = (.24756 - $MC_peaker) * mw_reductions_ces * 1000/12
gen welfare_loss_rectangle_sim_ces = (.24756 - $MC_peaker) * mw_reductions_sim_ces * 1000/12


foreach type in welfare_loss_rectangle welfare_loss_rectangle_sim welfare_loss_rectangle_ces welfare_loss_rectangle_sim_ces  welfare_loss_rectangle_non_event{
replace `type'  = 0 if hour < 14
replace `type' = 0 if hour >17
}



egen total_CS_losses = sum(welfare_losses+welfare_loss_rectangle), by(date)

//basically it takes all of the 5 minute increments within a day and adds them up. The assumption is that each event day has the same, which is fine

egen total_CS_losses_sim = sum(welfare_losses_sim+ welfare_loss_rectangle_sim), by(date)

egen total_welfare_losses_non_event = sum(welfare_loss_rectangle_non_event), by(date)

egen total_welfare_gains_1cent = sum(welfare_gains_1cent), by(date)
egen total_welfare_gains_1cent_sim = sum(welfare_gains_1cent_sim), by(date)


egen total_CS_losses_r1_tall_triangle = sum(welfare_losses_R1_tall_triangle + welfare_loss_rectangle), by(date) //uses the same rectangle as normal scenario b/c only 118 of reduction


///do for CES

egen total_CS_losses_ces = sum(welfare_losses_ces+welfare_loss_rectangle_ces), by(date)
egen total_CS_losses_sim_ces = sum(welfare_losses_sim_ces+ welfare_loss_rectangle_sim_ces), by(date)

egen total_welfare_gains_1cent_ces = sum(welfare_gains_1cent_ces), by(date)


//CALCULATE benefits from no longer needing to generate this power. Valued at peaker level

gen value = $MC_peaker * 1000 * mw_reductions /12
replace value = 0 if hour <14
replace value = 0 if hour > 17

gen value_simulated = $MC_peaker * 1000 * mw_reductions_sim /12
replace value_simulated = 0 if hour <14
replace value_simulated = 0 if hour > 17

gen value_ces = $MC_peaker * 1000 * mw_reductions_ces /12
replace value_ces = 0 if hour <14
replace value_ces = 0 if hour > 17

gen value_simulated_ces = $MC_peaker * 1000 * mw_reductions_sim_ces /12
replace value_simulated_ces = 0 if hour <14
replace value_simulated_ces = 0 if hour > 17



rename mw load

egen max_load = max(load), by(date)


*these 2 lines are for figure out a ratio and don't figure into final calculations - for a table
egen total_UTS_rectangle = sum( welfare_loss_rectangle_sim), by(date)
egen total_CS_lossestriangle= sum( welfare_losses), by(date)

keep date    event_day  total_CS_loss* mw_reductions mw_reductions_ces  mw_reductions_sim mw_reductions_sim_ces  total_welfare_gains_1cent total_welfare_gains_1cent_sim total_welfare_gains_1cent_ces    max_load  total_UTS_rectangle total_CS_lossestriangle total_CS_losses_r1_tall_triangle total_welfare_losses_non_event 
duplicates drop



gsort -max_load

save "stata data/max load for welfare calculations", replace
save "stata data/max load for welfare calculations $mw_reductions_value price_adder_$new_price_change_varlabel", replace



use "stata data/max load for welfare calculations", clear
egen summer_long_1cent_gains = sum(total_welfare_gains_1cent)
egen summer_long_1cent_gains_sim = sum(total_welfare_gains_1cent_sim)
egen summer_long_1cent_gains_sim_ces = sum(total_welfare_gains_1cent_ces)


gen dow = dow(date)
replace total_welfare_losses_non_event = 0 if inlist(dow,0,6)

egen summer_long_non_event_losses = sum(total_welfare_losses_non_event)


keep if event_day == 1



merge 1:1 date using "raw data/trigger and forcasted temps"
keep if _merge == 3

gsort -max_load

//generate 3 indicators for the 3 scenarios
gen trigger_day_keeper = 0
replace trigger_day_keeper = 1 if forecast_temps>=101

gen top_3 = 0
replace top_3 = 1 if _n<=3

gen all_days = 1




foreach day_scenario in all_days top_3 trigger_day_keeper  {

preserve
keep if `day_scenario' == 1

egen summer_CS_losses = sum(total_CS_losses)
egen summer_CS_losses_sim= sum(total_CS_losses_sim)

egen summer_CS_losses_ces = sum(total_CS_losses)
egen summer_CS_losses_sim_ces= sum(total_CS_losses_sim_ces)

*this is for R1's tall triangle question
egen summer_CS_losses_R1_tall_tri = sum(total_CS_losses_r1_tall_triangle)


keep mw_reductions* summer* summer_long_1cent_gains* summer_long_non_event_losses  
duplicates drop


foreach type in "" "_sim" "_sim_ces"{


gen per_mw_construction_cost`type' = 1159000*184/180 //from page 137 of cec report the 184/180 adjusts to 2016 data using PCCI from 

gen total_construction_costs`type' = mw_reductions`type'* per_mw_construction_cost`type'
gen construction_30_year`type' = total_construction_costs`type'/30

//o&M costs
//taken as $25.95 /kW/year from page 139 of CEC report

gen yearly_om_costs`type' = 25.95 * mw_reductions`type'*1000


//roughly this is mw reductions * hours (15*4) * emissions factor of peaker plant per mwh * conversion of pounds to short tons * social cost of carbon (50)
gen annual_value_co2`type' = mw_reductions`type' * 4*15 * 1239.3 * 0.000453592 * 50
*gen annual_value_co2`type' = mw_reductions`type' * 4*15 * 2180  * 0.000453592 * 50  //THis is if you want to do coal instead of NG - for R3


gen annual_value`type' =  yearly_om_costs`type' - summer_CS_losses`type'  + summer_long_1cent_gains`type' + annual_value_co2`type' // - summer_long_non_event_losses

***
*these are JUST to answer R1s question. It only works for default scenario (no sim) and it's the tall triangle
***
gen annual_value_tall_tri`type' = yearly_om_costs`type' - summer_CS_losses_R1_tall_tri  + summer_long_1cent_gains`type' + annual_value_co2`type'
***
*these are JUST to answer R1s question
***

//this high_cs is when I double the CS losses for a robustness check
gen annual_value`type'_highcs =  yearly_om_costs`type' - summer_CS_losses`type'*2  + summer_long_1cent_gains *2



gen required_cs_for_zero`type' =  yearly_om_costs`type'  + summer_long_1cent_gains`type'

gen separate_cs_benefits`type' =  yearly_om_costs`type'  + summer_long_1cent_gains`type' + annual_value_co2`type'
gen separate_cs_costs`type' =  - summer_CS_losses`type'

*I could separate the benefits from the costs and calclate separately? 


keep if _n == 1 //this is b/c I make the below calculations using 1 row, and the info needed is all duplicated
save "stata data//outcomes day_scenario `day_scenario' new price change $new_price_change_varlabel sensitivity scalar $sensitivity_scalar in progress.dta", replace

}
expand 30
gen discount = .03 //low discount rate chosens



pvvar annual_value discount, due(1)  //the due(1) makes it so the first period isn't discounted
pvvar annual_value_highcs discount, due(1)  //the due(1) makes it so the first period isn't discounted
pvvar annual_value_sim discount, due(1)  //the due(1) makes it so the first period isn't discounted
pvvar annual_value_sim_highcs discount, due(1)  //the due(1) makes it so the first period isn't discounted

pvvar annual_value_sim_ces discount, due(1)  //the due(1) makes it so the first period isn't discounted
pvvar annual_value_sim_ces_highcs discount, due(1)  //the due(1) makes it so the first period isn't discounted


pvvar required_cs_for_zero discount, due(1)  //the due(1) makes it so the first period isn't discounted
pvvar required_cs_for_zero_sim discount, due(1)  //the due(1) makes it so the first period isn't discounted
pvvar required_cs_for_zero_sim_ces discount, due(1)  //the due(1) makes it so the first period isn't discounted

pvvar separate_cs_benefits discount, due(1)  //the due(1) makes it so the first period isn't discounted
pvvar separate_cs_benefits_sim discount, due(1)  //the due(1) makes it so the first period isn't discounted
pvvar separate_cs_benefits_sim_ces discount, due(1)  //the due(1) makes it so the first period isn't discounted

pvvar separate_cs_costs discount, due(1)  //the due(1) makes it so the first period isn't discounted
pvvar separate_cs_costs_sim discount, due(1)  //the due(1) makes it so the first period isn't discounted
pvvar separate_cs_costs_sim_ces discount, due(1)  //the due(1) makes it so the first period isn't discounted


egen discounted_value = sum(annual_value_d1)
replace discounted_value = discounted_value  + total_construction_costs

egen discounted_value_highcs = sum(annual_value_highcs_d1)
replace discounted_value_highcs = discounted_value_highcs  + total_construction_costs

egen discounted_value_sim = sum(annual_value_sim_d1)
replace discounted_value_sim = discounted_value_sim  + total_construction_costs_sim

egen discounted_value_sim_highcs = sum(annual_value_sim_highcs_d1)
replace discounted_value_sim_highcs = discounted_value_sim_highcs  + total_construction_costs_sim

egen discounted_value_sim_ces = sum(annual_value_sim_ces_d1)
replace discounted_value_sim_ces = discounted_value_sim_ces  + total_construction_costs_sim_ces

egen discounted_value_sim_ces_highcs = sum(annual_value_sim_ces_highcs_d1)
replace discounted_value_sim_ces_highcs = discounted_value_sim_ces_highcs  + total_construction_costs_sim_ces


***
*these are JUST to answer R1s question
***
pvvar annual_value_tall_tri discount, due(1)  //the due(1) makes it so the first period isn't discounted
egen discounted_value_tall_tri = sum(annual_value_tall_tri_d1)
replace discounted_value_tall_tri = discounted_value_tall_tri + total_construction_costs





egen disval_separate_cs_benefits = sum(separate_cs_benefits_d1)

egen disval_separate_cs_benefits_sim = sum(separate_cs_benefits_sim_d1)
egen disval_separate_cs_bene_sim_ces = sum(separate_cs_benefits_sim_ces_d1) //note benefits shortened to bene

egen disval_separate_cs_costs = sum(separate_cs_costs_d1)

egen disval_separate_cs_costs_sim = sum(separate_cs_costs_sim_d1)
egen disval_separate_cs_costs_sim_ces = sum(separate_cs_costs_sim_ces_d1)


format discounted_value* disval_* annual* total* construction* yearly* %15.0fc
drop annual_value_*	annual_value_sim_* required_cs*

drop separate_cs_benefits_d1	separate_cs_benefits_sim_d1	separate_cs_costs_d1	separate_cs_costs_sim_d1 separate_cs_costs_sim_ces_d1 separate_cs_benefits_sim_ces_d1

duplicates drop

gen day_scenario = "`day_scenario'"

gen ratio = discounted_value_sim/discounted_value

gen simulated_price = $new_price_change

save "stata data//outcomes day_scenario `day_scenario' new price change $new_price_change_varlabel sensitivity scalar $sensitivity_scalar.dta", replace

restore
}
}
}






//a few other small things that are needed before the tables are made 


use if event_day == 1 using "stata data/A1 Standard regression dataset", clear

capture drop _merge


keep if total > 4000
drop if naics_most_2 == 22
drop if naics_most_2 == 51


drop if total > 50000

keep sa  naics_most_2 bin 
duplicates drop
save "stata data/saids used in main analysis", replace


use "stata data/saids used in main analysis", clear


merge 1:1 sa using "stata data\TOU_demographic 2"
keep if _merge == 3
drop _merge

drop if stopdate !=. 

drop naics2
merge m:1 sa using "stata data\TOU_demographic naics only"
drop if _merge == 2
drop _merge



keep sa rs  spid sp_min_date climate      bin naics_most_2
rename climate climate_str
encode climate_str, gen(climate)
drop climate_str

merge 1:1 sa using "stata data/PDP extract 2015_0713 full data"
drop if _merge == 2
gen pdp = 0
replace pdp = 1 if program_st == "Enrolled"


//I am defining defaulted as placed on to opt out rate. Based on document with andrew's info
//defaulted is anyone that has an action type. People who opt out are inactive status 

gen defaulted = 0
replace defaulted = 1 if action_type == "Default"
replace defaulted = 1 if action_type == "Affirm"
replace defaulted = 1 if action_type == "Optout"
replace defaulted = 1 if action_type == "Unenroll"

gen opt_out = 0
replace opt_out = 1 if action_type == "Optout"
replace opt_out = 1 if action_type == "Unenroll"
replace opt_out = . if _merge !=3
drop _merge

gen double sm_date = sp_min_date
format sm_date %td

capture drop sm_date_ym
gen sm_date_ym = mofd(sm_date)
format sm_date_ym %tm

drop if sm_date_ym > 635

save "stata data/A1 smart meter install date with interesting info", replace
